Generally, code blocks precede maps. To view code, click the “show” button.

Background and Significance

Transportation barriers pose significant challenges in providing consistent and timely care to patients. In a study using the National Health Interview Survey, nearly 2% of the United States population delayed medical care due to poor access to reliable transportation (Wolfe, McDonald, and Holmes 2020). Lack of reliable transportation can lead to missed appointments, delayed care, and inability to complete recommended treatment plans (Starbird et al. 2019; Syed, Gerber, and Sharp 2013).

Survivors of traumatic injury are a vulnerable patient population and face several barriers to accessing healthcare. Although the initial management of trauma is acute, ongoing survivorship of traumatic injuries require consistent follow up in the outpatient setting. Indeed, transportation barriers are cited in the literature as factors in influencing loss-to-follow up status (Jodoin et al. 2025). There is little insight in the literature regarding the spatial variation of transportation access and patient follow up. Understanding how transportation security varies across Chicago can give insight into how we may improve these barriers for our patient population.

Primary Means of Transportation

To understand the primary means of transportation around Chicago, data from the U.S. Census Bureau American Community Survey in 2023 was used for each census tracts aggregated over 5 years. Maps were constructed by clipping block groups into the Chicago city boundary. The census tracts and city boundary were in NAD 83 so they were reprojected to WGS 84 to make the geo-coded trauma locations and for the r5 routing engine. Block groups for the above analysis were also constructed similarly.

Population, race and ethnicity, transportation means to work, vehicle availability, and number of households below the poverty level were pulled from the 2023 ACS via the tigris package and converted to percentages.

Box plots were constructed to visualize the distribution of percentages of residents in each block group that do not have a car, drive to work, and use public transportation to get to work. While the percentage of people driving to work is fairly normal, the percentage who take public transit and the percentage who do not have a car is not. Therefore, I opted to display the data with Jenks breaks - this would allow fair comparison across graphs without having to standardize non-normal data.

I prefer the ggplot package because I feel the resolution to be better. ggplot does not have an automatic jenks function though, so i had to compute them manually with the ClassInt package in R and specifically extract the “breaks” from the classIntervals function. I then created a data frame with categorized data by cutting each column by the computed jenks breaks. I now had sorted data by jenks ready for ggplotting. While this took extra time and research, the map appears cleaner than tmap.

###--------------------------means of transportation-------------------------###
#quick exploration of means of transportation to work
ggplot(data = chi_tracts,
       mapping = aes(y = pdrive)) +
  geom_boxplot(coef = 1.5)

ggplot(data = chi_tracts,
       mapping = aes(y = ppubtrans)) +
  geom_boxplot(coef = 1.5)

ggplot(data = chi_tracts,
       mapping = aes(y = pwo_vehicel)) +
  geom_boxplot(coef = 1.5)

#percentage of driving to work is fairly normal but the percentage who take public transit
#and the percentage who do not have a car is not; will use jenks to display
### --- plot with ggplot --- ###
#compute jenks
jenks_breaks_pdrive <- classIntervals(chi_tracts$pdrive, n = 5, style = "jenks")$brks
jenks_breaks_ppubtrans <- classIntervals(chi_tracts$ppubtrans, n = 5, style = "jenks")$brks
jenks_breaks_pwo_vehicle <- classIntervals(chi_tracts$pwo_vehicel, n = 5, style = "jenks")$brks
#categorize data
chi_tracts_jenks <- chi_tracts %>%
  mutate(drive_cat = cut(pdrive, breaks = jenks_breaks_pdrive, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(pubtrans_cat = cut(ppubtrans, breaks = jenks_breaks_ppubtrans, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(wo_vehicle_cat = cut(pwo_vehicel, breaks = jenks_breaks_pwo_vehicle, include.lowest = TRUE, dig.lab = 5))
#plot with jenks
#percentage who drive
drive_map <- ggplot() +
  geom_sf(
    data = chi_tracts_jenks,
    aes(fill = drive_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Percentage") +
  labs(
    title = "Percentage Who Drive to Work",
    caption = "Source: ACS 2023 5-year estimates (Table B08301)") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.95, 0.80),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank())
#percentage who do not have a car
car_map <- ggplot() +
  geom_sf(
    data = chi_tracts_jenks,
    aes(fill = wo_vehicle_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Percentage") +
  labs(
    title = "Percentage Who Do Not Have a Car",
    caption = "Source: ACS 2023 5-year estimates (Table B25044)") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.95, 0.80),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank())

#percentage who do not have a car
pubtrans_map <- ggplot() +
  geom_sf(
    data = chi_tracts_jenks,
    aes(fill = pubtrans_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Percentage") +
  labs(
    title = "Percentage Who Take Public Transit to Work",
    caption = "Source: ACS 2023 5-year estimates (Table B08301)") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.95, 0.80),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank())

#interactive versions
drive_map_int <- tm_shape(chi_tracts_jenks) +
  tm_polygons(
    col = "drive_cat",
    palette = "OrRd",
    title = "Percentage",
    border.col = NA,
    lwd = 0,
    alpha = 0.8) +
  tm_layout(
    title = "Percentage Who Drive to Work",
    bg.color = "white",
    legend.position = c("right", "top")
  )
car_map_int <- tm_shape(chi_tracts_jenks) +
  tm_polygons(
    col = "wo_vehicle_cat",
    palette = "OrRd",
    title = "Percentage",
    border.col = NA,
    lwd = 0,
    alpha = 0.8) +
  tm_layout(
    title = "Percentage Who Do Not Have a Car",
    bg.color = "white",
    legend.position = c("right", "top")
  )

pubtrans_map_int <- tm_shape(chi_tracts_jenks) +
  tm_polygons(
    col = "pubtrans_cat",
    palette = "OrRd",
    title = "Percentage",
    border.col = NA,
    lwd = 0,
    alpha = 0.8) +
  tm_layout(
    title = "Percentage Who Take Public Transit to Work",
    bg.color = "white",
    legend.position = c("right", "top")
  )

Transportation Maps

Car Availability

Drives to Work

Public Transit to Work

Travel Time Analysis

Much of my travel time analysis was derived from the transportation chapter in Geocomputation with R (Lovelace, Nowosad, and Muenchow 2019).

Origins and Destinations

The unit of analysis chosen for travel time analysis is the census block group. Block groups have higher spatial resolution than census tracts, have fairly uniform populations (as compared to census tracts), but are more manageable for computational efforts than the highest spatial resolution - census blocks. The population centroid functions as the origin point for each block group. Population centroids approximate the geographic points of ’balance” in a block groups population density. Block groups with a population of 0 are excluded.

The destinations are level on trauma centers in Chicago as determined by the Illinois Department of Public Health. Level one trauma centers accept the highest acuity of trauma patients, and therefore will serve the population of interest - trauma patients with such sufficient injury that they need regular follow up with a trauma surgeon. I only included 9 trauma centers within or in close proximity of the city boundary by filtering the full, geo-coded trauma center list used in class. This included:

  • Advocate Christ Medical Center

  • Ascension Saint Francis

  • Mount Sinai Hospital

  • Advocate Illinois Masonic Medical Center

  • Advocate Lutheran General Hospital

  • John H. Stroger, Jr. Hospital of Cook County

  • Loyola University Medical Center

  • Northwestern Memorial Hospital

  • University of Chicago

###--------------------------origins and destinations------------------------###
#origins - centroids
USbgs_popcenter20 <- block_group2020
ILbgs_popcenter20 <- USbgs_popcenter20 %>%
  filter(STATEFP == 17) %>% # Illinois == 17
  mutate(GEOID = paste0(STATEFP, COUNTYFP, TRACTCE, BLKGRPCE)) # create a single tract GEOID similar to our other tract datasets
chi_bgs_popcenter20 <- ILbgs_popcenter20 %>%
  dplyr::select(GEOID, POPULATION, LONGITUDE, LATITUDE) %>%
  filter(GEOID %in% chi_bgs$GEOID) %>% # filter Illinois block groups on Chicago block groups only
  filter(POPULATION > 0) %>% 
  st_as_sf(coords = c("LONGITUDE", "LATITUDE"), # shapefile creation
           crs = 4326) # needs to be WGS 84
tmap_mode('view')
origins_map <- tm_shape(chi_bgs) +
  tm_borders() +
  tm_shape(chi_bgs_popcenter20) +
  tm_dots(fill = "blue",
          size = 0.25)
#destinations - trauma centers
trauma_geo$facility <- as.character(trauma_geo$facility)
adult_level_one <- trauma_geo %>% 
  filter(facility %in% c(
    'Advocate Christ Medical Center',
    'Ascension Saint Francis - Emergency Room',
    'Mount Sinai Hospital',
    'Advocate Illinois Masonic Medical Center',
    'Advocate Lutheran General Hospital',
    'John H. Stroger, Jr. Hospital of Cook County',
    'Loyola University Medical Center',
    'Northwestern Memorial Hospital',
    'University of Chicago Emergency Department and Trauma Center'))
destinations_map <- tm_shape(chiCity) +
  tm_borders() +
  tm_shape(adult_level_one) +
  tm_dots(fill = "darkgreen",
          size = 0.25)

Origins and Destinations Maps

Origins: Block Group Centroids

Destinations: Level One Trauma Centers

Routing Network Construction

The routing network was constructed from the “r5r” package in R (Pereira et al. 2021; Conway, Byrd, and Eggermond 2018; Conway, Byrd, and Van Der Linden 2017; Conway and Stewart 2019). This package uses the open source, “Rapid Realistic Routing on Real-world and Reimagined networks,” or R5, routing engine. This routing engine uses real-world street-netowrk data from OpenStreetMap and public transportation data in the form of General Transit Feed Specification (GTFS) paired with advance computational methods that account for uncertainty in public transit reliability to generate realistic estimations of public transit travel times.

Roadway and footpath data for Cook County were obtained from OpenStreetMap, and GTFS feeds for public transportation (timetables, routes, and service details) were downloaded from the CTA, Metra, and PACE websites. These inputs were used in r5r to build a multimodal routing network. Trip parameters included a maximum walking time of 30 minutes, maximum total travel time of 2 hours, and a departure window of November 19, 2025 at 09:00 ± 10 minutes. Origins were defined as population centroids of each census block group, and destinations as all level I trauma centers.

For each origin–destination pair, r5r generated multiple feasible public transit itineraries. Average travel time, wait time, walking time, number of transfers, and number of modal “legs” (walk, bus, train) were calculated by aggregating across all itineraries within the time window. Drive times were estimated separately using the road network, assuming typical traffic for the specified departure time.

Because this process produced a very large set of origin–destination combinations, each block group was ultimately assigned to a single trauma center based on the minimum drive time. The minimum drive time was chosen to approximate the trauma center that would be closest to each block group if resident was transported by EMS.

The following code depicts the above. This code is complex and computations took several hours given the magnitude of origin-destination pairs. Therefore, the code is displayed below for reference but will not run. The output is saved as an RData file for direct use in subsequent sections. A sample of the output of this file is shown below.

head(final_min_time)
###--------------------------------Data--------------------------------------###
# A road network data set from OpenStreetMap in .pbf format (mandatory)
# A public transport feed in GTFS.zip format
setwd("~/Desktop/GIS")
load('base_maps.RData')
load('origins-destinations.RData')

###--------------------------Route Construction------------------------------###
data_path <- "~/Desktop/GIS/transit_files"
#list.files(data_path)
r5r_network <- build_network(data_path = data_path)

###--------------------------------Test--------------------------------------###
#read dates
# read the GTFS feed
#gtfs_CTA <- read_gtfs('~/Desktop/GIS/transit_files/CTA.zip')
#gtfs_Metra <- read_gtfs('~/Desktop/GIS/transit_files/Metra.zip')
#gtfs_PACE <- read_gtfs('~/Desktop/GIS/transit_files/PACE.zip')

# get calendar dates (from calendar.txt)
#calendar_dates <- gtfs_CTA$calendar   

origin <- data.frame(id = '474 N Lakeshore',
                     lat = as.numeric(41.89108953920192),
                     lon = as.numeric(-87.61458686048464))

origin <- data.frame(id = 'Millenium Station',
                     lat = as.numeric(41.884589597182945),
                     lon = as.numeric(-87.6247049739245))

# set departure datetime input
mode <- c("WALK", 'TRANSIT')
max_walk_time <- 30 # minutes
max_trip_duration <- 120 # minutes
departure_datetime <- as.POSIXct("19-11-2025 09:00:00",
                                 format = "%d-%m-%Y %H:%M:%S")
ttm_test <- travel_time_matrix(
  r5r_network,
  origins = origin,
  destinations = ucm,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  max_trip_duration = max_trip_duration
)

# calculate detailed itineraries
ttm_test <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = ucm,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)

###-------------------------Transit Time Analysis----------------------------###
#create origins data frame
chi_bgcentroids$id <- chi_bgcentroids$GEOID
origin <- chi_bgcentroids %>% 
  select(id, geometry)
st_crs(origin) #needs to be WGS 84

#create destinations data frame
adult_level_one$id <- adult_level_one$facility
destination <- adult_level_one %>% 
  select(id, geometry)
st_crs(destination) #needs to be WGS 84

#UCM travel time matrix
ucm <- destination %>% 
  filter(id == 'University of Chicago Emergency Department and Trauma Center')
ttm_ucm <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = ucm,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Christ travel time matrix
christ <- destination %>% 
  filter(id == 'Advocate Christ Medical Center')
ttm_christ <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = christ,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Masonic travel time matrix
masonic <- destination %>% 
  filter(id == 'Advocate Illinois Masonic Medical Center')
ttm_masonic <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = masonic,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#County travel time matrix
county <- destination %>% 
  filter(id == 'John H. Stroger, Jr. Hospital of Cook County')
ttm_county <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = county,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#NWM travel time matrix
nwm <- destination %>% 
  filter(id == 'Northwestern Memorial Hospital')
ttm_nwm <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = nwm,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#St. Francis time matrix
sfh <- destination %>% 
  filter(id == 'Ascension Saint Francis - Emergency Room')
ttm_sfh <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = sfh,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Mount Sinai time matrix
msh <- destination %>% 
  filter(id == 'Mount Sinai Hospital')
ttm_msh <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = msh,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Loyola travel time matrix
loyola <- destination %>% 
  filter(id == 'Loyola University Medical Center')
ttm_loyola <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = loyola,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Lutheran travel time matrix
lutheran <- destination %>% 
  filter(id == 'Advocate Lutheran General Hospital')
ttm_lutheran <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = lutheran,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
###-------------------------Drive Time Analysis----------------------------###
#set mode to drive
mode <- c("CAR")
#UCM travel time matrix
dtm_ucm <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = ucm,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Christ travel time matrix
christ <- destination %>% 
  filter(id == 'Advocate Christ Medical Center')
dtm_christ <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = christ,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Masonic travel time matrix
masonic <- destination %>% 
  filter(id == 'Advocate Illinois Masonic Medical Center')
dtm_masonic <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = masonic,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#County travel time matrix
county <- destination %>% 
  filter(id == 'John H. Stroger, Jr. Hospital of Cook County')
dtm_county <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = county,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#NWM travel time matrix
nwm <- destination %>% 
  filter(id == 'Northwestern Memorial Hospital')
dtm_nwm <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = nwm,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#St. Francis time matrix
sfh <- destination %>% 
  filter(id == 'Ascension Saint Francis - Emergency Room')
dtm_sfh <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = sfh,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Mount Sinai time matrix
msh <- destination %>% 
  filter(id == 'Mount Sinai Hospital')
dtm_msh <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = msh,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Loyola travel time matrix
dtm_loyola <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = loyola,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Lutheran travel time matrix
dtm_lutheran <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = lutheran,
  mode = mode,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)

###--------------------Creating a Simplified Data Frame----------------------###
#below is a function that process the public transit data into one observation
#per origin-destination pair

process_ttm <- function(df) {
  
  df2 <- df %>%
    mutate(
      walk = if_else(mode == "WALK", 1, 0),
      transfer = if_else(wait > 0, 1, 0),
      bus = if_else(mode == "BUS", 1, 0),
      cta_train = if_else(mode == "SUBWAY", 1, 0),
      metra_train = if_else(mode == "RAIL", 1, 0)
    )
  
  walk_summary <- df2 %>%
    st_drop_geometry() %>%
    filter(mode == "WALK") %>%
    group_by(from_id, option) %>%
    summarise(
      total_walk_time  = sum(segment_duration, na.rm = TRUE),
      total_walk_distance = sum(distance, na.rm = TRUE),
      .groups = "drop"
    )
  
  wait_summary <- df2 %>%
    st_drop_geometry() %>%
    group_by(from_id, option) %>%
    summarise(
      wait_total = sum(wait, na.rm = TRUE),
      total_transfers = sum(transfer, na.rm = TRUE),
      total_bus_trips = sum(bus, na.rm = TRUE),
      total_cta_train_trips = sum(cta_train, na.rm = TRUE),
      total_metra_train_trips = sum(metra_train, na.rm = TRUE),
      number_of_legs = n(),
      .groups = "drop"
    )
  
  df3 <- df2 %>%
    left_join(walk_summary, by = c("from_id", "option")) %>%
    left_join(wait_summary, by = c("from_id", "option"))
  
  final <- df3 %>%
    st_drop_geometry() %>%
    group_by(from_id, to_id) %>%
    summarise(
      number_of_options = max(option),
      avg_tt = mean(total_duration, na.rm = TRUE),
      avg_td = mean(total_distance, na.rm = TRUE),
      avg_walk_time = mean(total_walk_time, na.rm = TRUE),
      avg_walk_dist = mean(total_walk_distance, na.rm = TRUE),
      avg_wait_time = mean(wait_total, na.rm = TRUE),
      avg_transfers = mean(total_transfers, na.rm = TRUE),
      avg_bus_trips = mean(total_bus_trips, na.rm = TRUE),
      avg_cta_train_trips = mean(total_cta_train_trips, na.rm = TRUE),
      avg_metra_train = mean(total_metra_train_trips, na.rm = TRUE),
      avg_legs = mean(number_of_legs, na.rm = TRUE),
      .groups = "drop"
    )
  final <- final %>% 
    rename(GEOID = from_id)
  
  return(final)
}
transit_ucm <- process_ttm(ttm_ucm)
transit_christ <- process_ttm(ttm_christ)
transit_county <- process_ttm(ttm_county)
transit_loyola <- process_ttm(ttm_loyola)
transit_lutheran <- process_ttm(ttm_lutheran)
transit_masonic <- process_ttm(ttm_masonic)
transit_msh <- process_ttm(ttm_msh)
transit_nwm <- process_ttm(ttm_nwm)
transit_sfh <- process_ttm(ttm_sfh)

#merged data
final_christ <- transit_christ %>% 
  full_join(dtm_christ_clean, by = "GEOID")
final_county <- transit_county %>% 
  full_join(dtm_county_clean, by = "GEOID")
final_loyola <- transit_loyola %>% 
  full_join(dtm_loyola_clean, by = "GEOID")
final_lutheran <- transit_lutheran %>% 
  full_join(dtm_lutheran_clean, by = "GEOID")
final_masonic <- transit_masonic %>% 
  full_join(dtm_masonic_clean, by = "GEOID")
final_msh <- transit_msh %>% 
  full_join(dtm_msh_clean, by = "GEOID")
final_nwm <- transit_nwm %>% 
  full_join(dtm_nwm_clean, by = "GEOID")
final_sfh <- transit_sfh %>% 
  full_join(dtm_sfh_clean, by = "GEOID")
final_ucm <- transit_ucm %>% 
  full_join(dtm_ucm_clean, by = "GEOID")

###------------------Creating closest drive time data frame--------------------###
#create mega data frame
final_combo <- bind_rows(final_christ,
          final_county,
          final_loyola,
          final_lutheran,
          final_masonic,
          final_msh,
          final_nwm,
          final_sfh,
          final_ucm)
#only keep closest by drive time
min_drive <- final_combo %>%
  group_by(GEOID) %>%
  slice_min(drive_time, with_ties = FALSE) %>%
  ungroup()
#combine with block group shapefiles
final_min_time <- chi_bgs %>% 
  left_join(min_drive, by = 'GEOID')

Travel Time Chloropleths

Block group level data were aggregated to the census tract level for more intuitive presentation. Travel time discrepancy between public transit and driving was calculated as the difference between the estimated transit time and drive time. Therefore, this could also be thought of as the “time saved by driving” rather than taking public transportation. A “weighted” minimum travel time was then derived by combining transit and drive times according to household car ownership: transit time was multiplied by the proportion of households without a vehicle, and drive time by the proportion with a vehicle. This approach assumes that residents with access to a car would use it to reach the nearest trauma center, while those without a car would rely on public transit. This may better approximate the average minimum travel time to the nearest trauma center among residents in a census tract.

Chloropleths are constructed in ggplot with jenks calculated as previously described

###----------------------------Time variables--------------------------------###
final_min_time <- final_min_time %>%
  mutate(
    travel_disc = avg_tt - drive_time,
    weighted_travel = avg_tt * (pwo_vehicel / 100) +
      drive_time * (1 - pwo_vehicel / 100), na.rm = T)
###-----------------Aggregate Travel Data to Tract Level---------------------###
cols_to_scale <- c('avg_tt', 'avg_td', 'avg_walk_time', 'avg_walk_dist',
                   'avg_wait_time', 'avg_transfers', 'avg_bus_trips', 'avg_cta_train_trips',
                   'avg_metra_train', 'avg_legs', 'drive_time', 'drive_distance',
                   'travel_disc', 'weighted_travel')
#scale up each average by block group population
tracts_travel <- final_min_time %>%
  mutate(across(all_of(cols_to_scale), ~ .x * total_popE)) 
#drop last charecter on GEOID to convert to tract GEOID
tracts_travel$GEOID <- substr(tracts_travel$GEOID, 1, nchar(tracts_travel$GEOID) - 1)
#summate scaled data
tracts_travel <- tracts_travel %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarise(
    total_popE = sum(total_popE, na.rm = TRUE),
    avg_tt = sum(avg_tt, na.rm = TRUE),
    avg_td = sum(avg_td, na.rm = TRUE),
    avg_walk_time = sum(avg_walk_time, na.rm = TRUE),
    avg_walk_dist = sum(avg_walk_dist, na.rm = TRUE),
    avg_wait_time = sum(avg_wait_time, na.rm = TRUE),
    avg_transfers = sum(avg_transfers, na.rm = TRUE),
    avg_bus_trips = sum(avg_bus_trips, na.rm = TRUE),
    avg_cta_train_trips = sum(avg_cta_train_trips, na.rm = TRUE),
    avg_metra_train = sum(avg_metra_train, na.rm = TRUE),
    avg_legs = sum(avg_legs, na.rm = TRUE),
    drive_time = sum(drive_time, na.rm = TRUE),
    drive_distance = sum(drive_distance, na.rm = TRUE),
    travel_disc = sum(travel_disc, na.rm = TRUE),
    weighted_travel = sum(weighted_travel, na.rm = TRUE),
    .groups = "drop")
#scale back down to average by dividing by total population
tracts_travel <- tracts_travel %>%
  mutate(across(all_of(cols_to_scale), ~ .x / total_popE)) 
#combine with shapefile
final_tracts <- chi_tracts %>% 
  left_join(tracts_travel, by = 'GEOID')
#compute jenks
jenks_breaks_transit <- classIntervals(final_tracts$avg_tt, n = 5, style = "jenks")$brks
jenks_breaks_drive <- classIntervals(final_tracts$drive_time, n = 5, style = "jenks")$brks
jenks_breaks_disc <- classIntervals(final_tracts$travel_disc, n = 5, style = "jenks")$brks
jenks_breaks_weight <- classIntervals(final_tracts$weighted_travel, n = 5, style = "jenks")$brks
#categorize data
travel_jenks <- final_tracts %>%
  mutate(transittime_cat = cut(avg_tt, breaks = jenks_breaks_transit, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(drivetime_cat = cut(drive_time, breaks = jenks_breaks_drive, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(disc_cat = cut(travel_disc, breaks = jenks_breaks_disc, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(weight_cat = cut(weighted_travel, breaks = jenks_breaks_weight, include.lowest = TRUE, dig.lab = 5)) 
#plot with jenks

Travel Time Maps

Public Transit

#transit time
transittime_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = transittime_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)") +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Estimated Public Transit Times",
    caption = "Aggregated Across All Multimodal Trips Within a 10-Minute Window at 09:00") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
transittime_map

Driving

#drive time
drivetime_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = drivetime_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)") +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Estimated Drive Times",
    caption = "Assuming Normal Traffic Patterns at 09:00") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
drivetime_map

Travel Discrepancy

#travel time discrepency time
disc_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = disc_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)") +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Time Saved by Driving",
    caption = "The difference between estimated public transit and driving times at 09:00") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
disc_map

Weighted Travel Time

weight_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = weight_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)") +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Weighted Travel Times",
    caption = "Weighted minimum travel times based on household car ownership") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
weight_map

Spatial Autocorrelation

Spatial clustering of extremes in drive times, public transit times, and weighted travel times was assessed via spatial autocorrelation. Global Moran’s I determined overall spatial correlation while local Moran’s I identified clusters of high and low spatial correlation. Briefly, Moran’s I statistic is the weighted average similarity in the measurement of interest between some observation and all it’s nearest neighbors divided by the sum of the sample variance. Nearest neighbors are determined by queen’s contiguity and measures of similarity via spatial lag with row-standardized weights (i.e. the average of the measurement of interest across all neighbors for a single observation).

Three census block groups with no population were excluded from prior analysis but will need to be included for spatial autocorrelation to avoid statistical artifact. Drive times and transit times for the centroids of each block group excluded were calculated and included in spatial autocorrelation measures. Given the populations of these block groups are zero, there is no data on household car ownership and therefore weighted travel times cannot be calculated. Their weighted travel times are approximated to be the mean weighted travel time.

Global Moran statistics for drive times, transit times, and weighted travel times all indicated spatial autocorrelation (pseudo-p < 0.01). Local Moran’s I identified clustered areas of high access (travel times less than one standard deviation from the mean and a Moran’s I > 0) and areas of low access (travel times greater than one standard deviation from the mean and a Moran’s I > 0). Scatter plots and maps below. Note that there was one census tract with significant dispersion in the weighted travel times where the tract experiences low access but is surrounded by high access.

#some tracts have no drive time because they were excluded from the routing network
#as they had no population - O'Hare, Midway, and a small Railroad tract
#these tracts are replaced with manually computed travel times
spatial_final_tracts <- final_tracts
#O'Hare
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980100, "drive_time"] <- 16.67
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980100, "avg_tt"] <- 54.77
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980100, "weighted_travel"] <- mean(
  spatial_final_tracts$weighted_travel, na.rm = T)
#Railway Tract
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031381700, "drive_time"] <- 10
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031381700, "avg_tt"] <- 30.55
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031381700, "weighted_travel"] <- mean(
  spatial_final_tracts$weighted_travel, na.rm = T)
#Midway
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980000, "drive_time"] <- 14.33
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980000, "avg_tt"] <- 37
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980000, "weighted_travel"] <- mean(
  spatial_final_tracts$weighted_travel, na.rm = T)
###-------------------------Spatial Neighbors--------------------------------###
#queen contiguity
chi_tracts_q <- poly2nb(chi_tracts, queen = TRUE) #average 6.55 links
##spatial weights matrix
#row-standardized weights for spatial neighbors
tractweights_q <- nb2listw(chi_tracts_q, style = "W")

Moran’s I and Local Moran’s I Calculations

Drive Time

#spatial lag for drive time
spatial_final_tracts$lag_drive_q <- lag.listw(tractweights_q, spatial_final_tracts$drive_time)
#plot spatial lag for drive time
ggplot(spatial_final_tracts, aes(x = drive_time, y = lag_drive_q)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Estimated Drive Time (min) to nearest trauma center by census tract, Chicago",
       x = "Drive Time (min)",
       y = "Spatial lag, Drive Time (min)", 
       caption = "Spatial relationships based on queens-case polygon contiguity.")

### --- standardized data --- ###
mean_drivetime = mean(spatial_final_tracts$drive_time)
sd_drivetime = sqrt(var(spatial_final_tracts$drive_time))
mean_lag_drive_q = mean(spatial_final_tracts$lag_drive_q)
sd_lag_drive_q = sqrt(var(spatial_final_tracts$lag_drive_q))
spatial_final_tracts$std_drivetime <- (spatial_final_tracts$drive_time - mean_drivetime)/sd_drivetime
spatial_final_tracts$std_lag_drive_q <- (spatial_final_tracts$lag_drive_q - mean_lag_drive_q)/sd_lag_drive_q
#map of high stadard deviation census tracts
tm_shape(spatial_final_tracts) +
  tm_borders(col = "black") +
  tm_shape(spatial_final_tracts[(spatial_final_tracts$std_drivetime > 2 & spatial_final_tracts$std_lag_drive_q > 1),]) +
  tm_polygons(col = "red",
              border.col = "black") +
  tm_layout(main.title = "Drive Time High Std Deviation",
            main.title.size = 0.7,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
### --- Moran's I --- ###
set.seed(1234)
queen_rwstd_drive <- moran.mc(spatial_final_tracts$drive_time,
                             listw = tractweights_q, 
                             nsim = 999) 
queen_rwstd_drive
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  spatial_final_tracts$drive_time 
## weights: tractweights_q  
## number of simulations + 1: 1000 
## 
## statistic = 0.9184, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
moran.plot(spatial_final_tracts$drive_time, listw=tractweights_q)

### --- Local Moran's I --- ###
set.seed(1234)
local_drive <- localmoran_perm(
  spatial_final_tracts$std_drivetime, 
  tractweights_q, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
#data frame with local I
local_drive_df <- spatial_final_tracts %>%
  select(GEOID, std_drivetime,std_lag_drive_q) %>%
  bind_cols(local_drive)
#clusters data frame
drive_clusters <- local_drive_df %>%
  mutate(local_cluster = case_when(
    p_i >= 0.01 ~ "Not significant",
    std_drivetime > 0 & local_i > 0 ~ "Low access",
    std_drivetime > 0 & local_i < 0 ~ "Low-high",
    std_drivetime < 0 & local_i > 0 ~ "High access",
    std_drivetime < 0 & local_i < 0 ~ "High-low"))
#plot local Moran's I
color_values <- c(`Low access` = "red", 
                  `Low-high` = "pink", 
                  `High access` = "blue", 
                  `High-low` = "lightblue", 
                  `Not significant` = "white")

ggplot(drive_clusters, aes(x = std_drivetime, 
                                      y = std_lag_drive_q,
                                      fill = local_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Drive time (z-score)",
       y = "Spatial lag of drive time (z-score)",
       fill = "Cluster type")

#mapping access by car
drive_access <- tm_shape(drive_clusters) +
  tm_polygons(col = "local_cluster",
              palette = color_values,
              title = "Level of Access",
              textNA = "Undefined") +
  tm_layout(main.title = "Drive Times to Nearest Level One Trauma Center, Chicago (2023)",
            main.title.size = 0.8,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values'), 'textNA'
## (rename to 'label.na') to fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

Transit Time

###-------------------------Transit time Analysis------------------------------###

#spatial lag for transit time
spatial_final_tracts$lag_transit_q <- lag.listw(tractweights_q, spatial_final_tracts$avg_tt)
#plot spatial lag for transit time
ggplot(spatial_final_tracts, aes(x = avg_tt, y = lag_transit_q)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Estimated Transit Time (min) to nearest trauma center by census tract, Chicago",
       x = "Transit Time (min)",
       y = "Spatial lag, Transit Time (min)", 
       caption = "Spatial relationships based on queens-case polygon contiguity.")

### --- standardized data --- ###
mean_transit = mean(spatial_final_tracts$avg_tt)
sd_transit = sqrt(var(spatial_final_tracts$avg_tt))
mean_lag_transit_q = mean(spatial_final_tracts$lag_transit_q)
sd_lag_transit_q = sqrt(var(spatial_final_tracts$lag_transit_q))
spatial_final_tracts$std_transit <- (spatial_final_tracts$avg_tt - mean_transit)/sd_transit
spatial_final_tracts$std_lag_transit_q <- (spatial_final_tracts$lag_transit_q - mean_lag_transit_q)/sd_lag_transit_q
#map of high stadard deviation census tracts
tm_shape(spatial_final_tracts) +
  tm_borders(col = "black") +
  tm_shape(spatial_final_tracts[(spatial_final_tracts$std_transit > 2 & spatial_final_tracts$std_lag_transit_q > 1),]) +
  tm_polygons(col = "red",
              border.col = "black") +
  tm_layout(main.title = "Transit Time High Std Deviation",
            main.title.size = 0.7,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
### --- Moran's I --- ###
set.seed(1235)
queen_rwstd_transit <- moran.mc(spatial_final_tracts$avg_tt,
                              listw = tractweights_q, 
                              nsim = 999) 
queen_rwstd_transit
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  spatial_final_tracts$avg_tt 
## weights: tractweights_q  
## number of simulations + 1: 1000 
## 
## statistic = 0.89724, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
moran.plot(spatial_final_tracts$avg_tt, listw=tractweights_q)

### --- Local Moran's I --- ###
set.seed(1235)
local_transit <- localmoran_perm(
  spatial_final_tracts$std_transit, 
  tractweights_q, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
#data frame with local I
local_transit_df <- spatial_final_tracts %>%
  select(GEOID, std_transit,std_lag_transit_q) %>%
  bind_cols(local_transit)
#clusters data frame
transit_clusters <- local_transit_df %>%
  mutate(local_cluster = case_when(
    p_i >= 0.01 ~ "Not significant",
    std_transit > 0 & local_i > 0 ~ "Low access",
    std_transit > 0 & local_i < 0 ~ "Low-high",
    std_transit < 0 & local_i > 0 ~ "High access",
    std_transit < 0 & local_i < 0 ~ "High-low"))
#plot local Moran's I
ggplot(transit_clusters, aes(x = std_transit, 
                           y = std_lag_transit_q,
                           fill = local_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Transit Time (z-score)",
       y = "Spatial lag of transit time (z-score)",
       fill = "Cluster type")

#mapping access by public transit
transit_access <- tm_shape(transit_clusters) +
  tm_polygons(col = "local_cluster",
              palette = color_values,
              title = "Level of Access",
              textNA = "Undefined") +
  tm_layout(main.title = "Transit Travel Times to Nearest Level One Trauma Center, Chicago (2023)",
            main.title.size = 0.8,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values'), 'textNA'
## (rename to 'label.na') to fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

Weighted Travel Time

###----------------------Weighted Travel Time Analysis-----------------------###

#spatial lag for wieghted travel
spatial_final_tracts$lag_weight_q <- lag.listw(tractweights_q, spatial_final_tracts$weighted_travel)
#plot spatial lag for weighted time
ggplot(spatial_final_tracts, aes(x = weighted_travel, y = lag_weight_q)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Weighted Travel Time (min) to nearest trauma center by census tract, Chicago",
       x = "Weighted Travel Time (min)",
       y = "Spatial lag, Weighted Travel Time (min)", 
       caption = "Spatial relationships based on queens-case polygon contiguity.")

### --- standardized data --- ###
mean_weight = mean(spatial_final_tracts$weighted_travel)
sd_weight = sqrt(var(spatial_final_tracts$weighted_travel))
mean_lag_weight_q = mean(spatial_final_tracts$lag_weight_q)
sd_lag_weight_q = sqrt(var(spatial_final_tracts$lag_weight_q))
spatial_final_tracts$std_weight <- (spatial_final_tracts$weighted_travel - mean_weight)/sd_weight
spatial_final_tracts$std_lag_weight_q <- (spatial_final_tracts$lag_weight_q - mean_lag_weight_q)/sd_lag_weight_q
#map of high stadard deviation census tracts
tm_shape(spatial_final_tracts) +
  tm_borders(col = "black") +
  tm_shape(spatial_final_tracts[(spatial_final_tracts$std_weight > 2 & spatial_final_tracts$std_lag_weight_q > 1),]) +
  tm_polygons(col = "red",
              border.col = "black") +
  tm_layout(main.title = "Weighted Travel Time High Std Deviation",
            main.title.size = 0.7,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
### --- Moran's I --- ###
set.seed(1236)
queen_rwstd_weight <- moran.mc(spatial_final_tracts$weighted_travel,
                                listw = tractweights_q, 
                                nsim = 999) 
queen_rwstd_weight
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  spatial_final_tracts$weighted_travel 
## weights: tractweights_q  
## number of simulations + 1: 1000 
## 
## statistic = 0.79436, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
moran.plot(spatial_final_tracts$weighted_travel, listw=tractweights_q)

### --- Local Moran's I --- ###
set.seed(1236)
local_weight <- localmoran_perm(
  spatial_final_tracts$std_weight, 
  tractweights_q, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
#data frame with local I
local_weight_df <- spatial_final_tracts %>%
  select(GEOID, std_weight,std_lag_weight_q) %>%
  bind_cols(local_weight)
#clusters data frame
weight_clusters <- local_weight_df %>%
  mutate(local_cluster = case_when(
    p_i >= 0.01 ~ "Not significant",
    std_weight > 0 & local_i > 0 ~ "Low access",
    std_weight > 0 & local_i < 0 ~ "Low-high",
    std_weight < 0 & local_i > 0 ~ "High access",
    std_weight < 0 & local_i < 0 ~ "High-low"))
#plot local Moran's I
ggplot(weight_clusters, aes(x = std_weight, 
                             y = std_lag_weight_q,
                             fill = local_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Weighted Travel Time (z-score)",
       y = "Spatial lag of weighted travel time (z-score)",
       fill = "Cluster type")

#mapping access by car
weight_access <- tm_shape(weight_clusters) +
  tm_polygons(col = "local_cluster",
              palette = color_values,
              title = "Level of Access",
              textNA = "Undefined") +
  tm_layout(main.title = "Weighted Travel Times to Nearest Level One Trauma Center, Chicago (2023)",
            main.title.size = 0.8,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'palette' (rename to 'values'), 'textNA'
## (rename to 'label.na') to fill.scale = tm_scale(<HERE>).
## ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
##   a list: 'fill.scale = list(<scale1>, <scale2>, ...)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`

Mapping Access

Drive Time

drive_access

Transit Time

transit_access

Weighted Travel Time

weight_access

Demographics and Travel

atlas_data$GEOID <- as.character(atlas_data$GEOID)
final_tracts <- final_tracts %>% 
  left_join(atlas_data, by = 'GEOID')
final_tracts$transit_access_level <- transit_clusters$local_cluster
final_tracts$weight_access_level <- weight_clusters$local_cluster
tracts_access <- final_tracts %>%
  filter(transit_access_level != "Not significant") %>%
  filter(total_popE.x > 0) %>%
  mutate(psenior_pov = ifelse(is.infinite(psenior_pov), 0, psenior_pov)) %>% 
  mutate(psenior_disab = ifelse(is.infinite(psenior_disab), 0, psenior_disab)) %>% 
  mutate(psenior_povdisab = ifelse(is.infinite(psenior_povdisab), 0, psenior_povdisab))
tracts_access_transit <- tracts_access %>%
  st_drop_geometry() %>% 
  mutate(transit_access_level = if_else(
    transit_access_level == 'Low access', 1, 0),
    transit_access_level = factor(
      transit_access_level,
      levels = c(0,1),
      labels = c("High Access", "Low Access")
    ))
tracts_access_weight <- tracts_access %>%
  st_drop_geometry() %>% 
  mutate(weight_access_level = if_else(
    weight_access_level == 'Low access', 1, 0),
    weight_access_level = factor(
      weight_access_level,
      levels = c(0,1),
      labels = c("High Access", "Low Access")
    ))

Demographic and Transit-Time Characteristics of Census Tracts

The below tables represent the differences between high and low access groups across several demographic measures.

Transit Times Demographics

demographics_table_transit <-
  tracts_access_transit %>%
  select(transit_access_level,
         total_popE.x, pwhite, pblack, pasian, phisp,
         psenior, psenior_pov, psenior_disab,
         psenior_povdisab, ppov) %>%
  tbl_summary(
    by = transit_access_level,
    missing = "no",
    statistic = list(
      everything() ~ "{mean} ({sd})",
      total_popE.x ~ "{sum}"
    ),
    label = list(
      total_popE.x = "Total Population",
      pwhite = "Non-Hispanic White",
      pblack = "Non-Hispanic Black",
      pasian = "Non-Hispanic Asian",
      phisp = "Hispanic",
      psenior = "Age 65+",
      psenior_pov = "65+ and Under Poverty Level",
      psenior_disab = "65+ and Disabled",
      psenior_povdisab = "65+ Under Poverty Level and Disabled",
      ppov = "Under Poverty Level"
    )
  ) %>%
  add_p() %>%
  as_gt() %>%
  tab_row_group(
    group = md("**Socioeconomic Indicators (%)**"),
    rows = variable %in% c("ppov")) %>% 
  tab_row_group(
    group = md("**Senior Characteristics (%)**"),
    rows = variable %in% c("psenior", "psenior_pov", 
                           "psenior_disab", "psenior_povdisab")) %>% 
  tab_row_group(
    group = md("**Demographics (%)**"),
    rows = variable %in% c("total_popE.x", 
                           "pwhite", "pblack", "pasian", "phisp")
  )
demographics_table_transit
Characteristic High Access
N = 126
1
Low Access
N = 116
1
p-value2
Demographics (%)
Total Population 410,416 488,626 <0.001
Non-Hispanic White 47 (30) 32 (27) <0.001
Non-Hispanic Black 27 (33) 22 (34) <0.001
Non-Hispanic Asian 8 (8) 4 (6) <0.001
Hispanic 14 (20) 39 (28) <0.001
Senior Characteristics (%)
Age 65+ 6.2 (5.3) 9.4 (3.7) <0.001
65+ and Under Poverty Level 21 (39) 19 (28) 0.023
65+ and Disabled 61 (78) 59 (71) 0.055
65+ Under Poverty Level and Disabled 32 (80) 19 (39) 0.009
Socioeconomic Indicators (%)
Under Poverty Level 11 (13) 11 (10) 0.025
1 Sum; Mean (SD)
2 Wilcoxon rank sum test

Weighted Travel Times Demographics

demographics_table_weight <-
  tracts_access_weight %>%
  select(weight_access_level,
         total_popE.x, pwhite, pblack, pasian, phisp,
         psenior, psenior_pov, psenior_disab,
         psenior_povdisab, ppov) %>%
  tbl_summary(
    by = weight_access_level,
    missing = "no",
    statistic = list(
      everything() ~ "{mean} ({sd})",
      total_popE.x ~ "{sum}"
    ),
    label = list(
      total_popE.x = "Total Population",
      pwhite = "Non-Hispanic White",
      pblack = "Non-Hispanic Black",
      pasian = "Non-Hispanic Asian",
      phisp = "Hispanic",
      psenior = "Age 65+",
      psenior_pov = "65+ and Under Poverty Level",
      psenior_disab = "65+ and Disabled",
      psenior_povdisab = "65+ Under Poverty Level and Disabled",
      ppov = "Under Poverty Level"
    )
  ) %>%
  add_p() %>%
  as_gt() %>%
  tab_row_group(
    group = md("**Socioeconomic Indicators (%)**"),
    rows = variable %in% c("ppov")) %>% 
  tab_row_group(
    group = md("**Senior Characteristics (%)**"),
    rows = variable %in% c("psenior", "psenior_pov", 
                           "psenior_disab", "psenior_povdisab")) %>% 
  tab_row_group(
    group = md("**Demographics (%)**"),
    rows = variable %in% c("total_popE.x", 
                           "pwhite", "pblack", "pasian", "phisp")
  )
demographics_table_weight
Characteristic High Access
N = 170
1
Low Access
N = 72
1
p-value2
Demographics (%)
Total Population 604,426 294,616 0.002
Non-Hispanic White 47 (29) 23 (23) <0.001
Non-Hispanic Black 22 (31) 31 (39) 0.5
Non-Hispanic Asian 8 (8) 3 (4) <0.001
Hispanic 20 (23) 41 (30) <0.001
Senior Characteristics (%)
Age 65+ 7.2 (5.3) 8.9 (3.5) <0.001
65+ and Under Poverty Level 19 (35) 22 (33) 0.009
65+ and Disabled 58 (69) 66 (88) 0.14
65+ Under Poverty Level and Disabled 26 (69) 25 (49) <0.001
Socioeconomic Indicators (%)
Under Poverty Level 10 (12) 14 (11) <0.001
1 Sum; Mean (SD)
2 Wilcoxon rank sum test

Transit Time Travel Metrics

transit_table_unweighted <-
  tracts_access_transit %>%
  select(transit_access_level,
         total_popE.x,
         drive_time, avg_tt,
         ppubtrans, pdrive, pwo_vehicel,
         avg_walk_time, avg_wait_time,
         avg_transfers, avg_legs,
         avg_bus_trips, avg_cta_train_trips,
         EKW_2024, RITB_2022) %>%
  tbl_summary(
    by = transit_access_level,
    missing = "no",
    statistic = list(
      everything() ~ "{mean} ({sd})",
      total_popE.x ~ "{sum}"
    ),
    label = list(
      total_popE.x = "Total Population",
      drive_time = "Drive Time",
      avg_tt = "Transit Time",
      ppubtrans = "Take Public Transit to Work",
      pdrive  = "Drive Alone to Work",
      pwo_vehicel = "No Vehicle Household",
      avg_walk_time = "Walk Time",
      avg_wait_time = "Wait Time",
      avg_transfers = "Transfers",
      avg_legs = "Journey Legs",
      avg_bus_trips = "Bus Trips",
      avg_cta_train_trips = "CTA Train Trips",
      EKW_2024 = "Walkability Access",
      RITB_2022 = "Transportation Burden"
    )
  ) %>%
  add_p() %>%
  as_gt() %>%
  tab_row_group(
    group = md("**Neighbourhood Characteristics (Percentile)**"),
    rows = variable %in% c("EKW_2024", "RITB_2022")) %>%
  tab_row_group(
    group = md("**Transit Trip Structure (Average Count per Trip)**"),
    rows = variable %in% c("avg_transfers",
                           "avg_legs",
                           "avg_bus_trips",
                           "avg_cta_train_trips")) %>% 
  tab_row_group(
    group = md("**Access Mode (%)**"),
    rows = variable %in% c("ppubtrans", "pdrive", "pwo_vehicel")) %>%
  tab_row_group(
    group = md("**Transit Time Components (minutes)**"),
    rows = variable %in% c("drive_time", "avg_tt",
                           "avg_walk_time", "avg_wait_time"))

transit_table_unweighted
Characteristic High Access
N = 126
1
Low Access
N = 116
1
p-value2
Transit Time Components (minutes)
Drive Time 5 (2) 21 (4) <0.001
Transit Time 19 (4) 64 (8) <0.001
Walk Time 12.1 (3.2) 17.8 (6.5) <0.001
Wait Time 2.1 (1.1) 8.9 (4.1) <0.001
Access Mode (%)
Take Public Transit to Work 25 (10) 12 (7) <0.001
Drive Alone to Work 33 (17) 63 (9) <0.001
No Vehicle Household 36 (14) 14 (11) <0.001
Transit Trip Structure (Average Count per Trip)
Transfers 1.14 (0.32) 2.31 (0.53) <0.001
Journey Legs 3.28 (0.63) 5.56 (1.02) <0.001
Bus Trips 1.00 (0.38) 2.07 (0.55) <0.001
CTA Train Trips 0.14 (0.25) 0.20 (0.30) 0.010
Neighbourhood Characteristics (Percentile)
Walkability Access 15.31 (1.77) 14.47 (1.40) <0.001
Transportation Burden 45 (21) 57 (25) <0.001
Total Population 410,416 488,626 <0.001
1 Sum; Mean (SD)
2 Wilcoxon rank sum test

Weighted Travel Times Travel Metrics

transit_table_weighted <-
  tracts_access_weight %>%
  select(weight_access_level,
         total_popE.x,
         drive_time, avg_tt,
         ppubtrans, pdrive, pwo_vehicel,
         avg_walk_time, avg_wait_time,
         avg_transfers, avg_legs,
         avg_bus_trips, avg_cta_train_trips,
         EKW_2024, RITB_2022) %>%
  tbl_summary(
    by = weight_access_level,
    missing = "no",
    statistic = list(
      everything() ~ "{mean} ({sd})",
      total_popE.x ~ "{sum}"
    ),
    label = list(
      total_popE.x = "Total Population",
      drive_time = "Drive Time",
      avg_tt = "Transit Time",
      ppubtrans = "Take Public Transit to Work",
      pdrive  = "Drive Alone to Work",
      pwo_vehicel = "No Vehicle Household",
      avg_walk_time = "Walk Time",
      avg_wait_time = "Wait Time",
      avg_transfers = "Transfers",
      avg_legs = "Journey Legs",
      avg_bus_trips = "Bus Trips",
      avg_cta_train_trips = "CTA Train Trips",
      EKW_2024 = "Walkability Access",
      RITB_2022 = "Transportation Burden"
    )
  ) %>%
  add_p() %>%
  as_gt() %>%
  tab_row_group(
    group = md("**Neighbourhood Characteristics (Percentile)**"),
    rows = variable %in% c("EKW_2024", "RITB_2022")) %>%
  tab_row_group(
    group = md("**Transit Trip Structure (Average Count per Trip)**"),
    rows = variable %in% c("avg_transfers",
                           "avg_legs",
                           "avg_bus_trips",
                           "avg_cta_train_trips")) %>% 
  tab_row_group(
    group = md("**Access Mode (%)**"),
    rows = variable %in% c("ppubtrans", "pdrive", "pwo_vehicel")) %>%
  tab_row_group(
    group = md("**Transit Time Components (minutes)**"),
    rows = variable %in% c("drive_time", "avg_tt",
                           "avg_walk_time", "avg_wait_time"))

transit_table_weighted
Characteristic High Access
N = 170
1
Low Access
N = 72
1
p-value2
Transit Time Components (minutes)
Drive Time 9 (6) 23 (4) <0.001
Transit Time 31 (21) 65 (8) <0.001
Walk Time 14.4 (6.1) 15.8 (4.9) 0.004
Wait Time 4.3 (4.5) 7.9 (3.3) <0.001
Access Mode (%)
Take Public Transit to Work 21 (11) 13 (7) <0.001
Drive Alone to Work 41 (21) 62 (9) <0.001
No Vehicle Household 29 (17) 16 (12) <0.001
Transit Trip Structure (Average Count per Trip)
Transfers 1.48 (0.69) 2.24 (0.51) <0.001
Journey Legs 3.94 (1.35) 5.40 (0.96) <0.001
Bus Trips 1.31 (0.68) 1.99 (0.51) <0.001
CTA Train Trips 0.17 (0.28) 0.18 (0.28) 0.3
Neighbourhood Characteristics (Percentile)
Walkability Access 15.18 (1.66) 14.25 (1.45) <0.001
Transportation Burden 44 (22) 66 (22) <0.001
Total Population 604,426 294,616 0.002
1 Sum; Mean (SD)
2 Wilcoxon rank sum test

Implications for Local Patients

Turning attention to the University of Chicago Patient Population

transit_ucm <- transit_ucm %>% 
  left_join(dtm_ucm_clean, by = "GEOID")
transit_ucm <- chi_bgs %>%
  left_join(transit_ucm, by = "GEOID")
ucm <- adult_level_one %>% 
  filter(facility == "University of Chicago Emergency Department and Trauma Center")
travel_breaks_ucm <- transit_ucm %>%
  mutate(transittime_cat = cut(
      avg_tt,
      breaks = c(0, 30, 60, Inf),
      include.lowest = TRUE,
      right = TRUE,
      labels = c("< 30 minutes", "30–60 minutes", "> 60 minutes"))) %>% 
  mutate(benefit = case_when(
    avg_tt <= 30 ~ 0,
    avg_tt >= 30 & drive_time <= 30 ~1,
    avg_tt >= 30 & drive_time >= 30 ~0),
    benefit = factor(
      benefit,
      levels = c(0,1),
      labels = c("Likely Low Impact", "Likely High Impact")
    ))

Potential Benefit of a Rideshare Program

Transit time to UChicago

ucm_transittime <- ggplot() +
  geom_sf(
    data = travel_breaks_ucm,
    aes(fill = transittime_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)") +
  geom_sf(
    data = ucm,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Estimated Public Transit Times",
    caption = "Aggregated Across All Multimodal Trips Within a 10-Minute Window at 09:00") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
ucm_transittime

Potential Benefit

ucm_benefit <- ggplot() +
  geom_sf(
    data = travel_breaks_ucm,
    aes(fill = benefit),
    color = NA) +
  scale_fill_manual(
    values = c(
      "Likely Low Impact" = "lightblue",
      "Likely High Impact" = "navyblue"),
    name = "Benefit") +
  geom_sf(
    data = ucm,
    color = "navyblue",
    fill = "navyblue",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Potential Benefit of Rideshare Program",
    caption = "Determined by 30-minute Time Cut Off") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
ucm_benefit

References

Conway, Matthew Wigginton, Andrew Byrd, and Michael van Eggermond. 2018. “Accounting for Uncertainty and Variation in Accessibility Metrics for Public Transport Sketch Planning.” Journal of Transport and Land Use 11 (1). https://doi.org/10.5198/jtlu.2018.1074.
Conway, Matthew Wigginton, Andrew Byrd, and Marco Van Der Linden. 2017. “Evidence-Based Transit and Land Use Sketch Planning Using Interactive Accessibility Methods on Combined Schedule and Headway-Based Networks.” Transportation Research Record: Journal of the Transportation Research Board 2653 (1): 45–53. https://doi.org/10.3141/2653-06.
Conway, Matthew Wigginton, and Anson F. Stewart. 2019. “Getting Charlie Off the MTA: A Multiobjective Optimization Method to Account for Cost Constraints in Public Transit Accessibility Metrics.” International Journal of Geographical Information Science 33 (9): 1759–87. https://doi.org/10.1080/13658816.2019.1605075.
Jodoin, Zachary, Daanish Sheikh, Cameron Atkinson, Loc Uyen Vo, Alvaro Moreira, Christina Brady, and Boris Zelle. 2025. “Paediatric ballistic fracture patients: who has poor follow-up and why?” International Orthopaedics 49 (6): 1451–60. https://doi.org/10.1007/s00264-025-06506-3.
Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2019. Chapter 13 Transportation | Geocomputation with R. https://r.geocompx.org/transport.html.
Pereira, Rafael H. M., Marcus Saraiva, Daniel Herszenhut, Carlos Kaue Vieira Braga, and Matthew Wigginton Conway. 2021. “R5r: Rapid Realistic Routing on Multimodal Transport Networks with R5 in R.” Findings, March. https://doi.org/10.32866/001c.21262.
Starbird, Laura E., Caitlin DiMaina, Chun-An Sun, and Hae-Ra Han. 2019. “A Systematic Review of Interventions to Minimize Transportation Barriers Among People with Chronic Diseases.” Journal of Community Health 44 (2): 400–411. https://doi.org/10.1007/s10900-018-0572-3.
Syed, Samina T., Ben S. Gerber, and Lisa K. Sharp. 2013. “Traveling Towards Disease: Transportation Barriers to Health Care Access.” Journal of Community Health 38 (5): 976–93. https://doi.org/10.1007/s10900-013-9681-1.
Wolfe, Mary K., Noreen C. McDonald, and G. Mark Holmes. 2020. “Transportation Barriers to Health Care in the United States: Findings from the National Health Interview Survey, 19972017.” American Journal of Public Health 110 (6): 815–22. https://doi.org/10.2105/AJPH.2020.305579.